TP 1: Regresión lineal

Diagnóstico y Evaluación de Modelos de Regresión Lineal

Dataset

Los datos con los que se trabajará en este TP provienen de la 3° Encuesta Mundial de Salud Escolar (EMSE) provistos por el Ministerio de Salud link de la República Argentina. Esta encuesta trata sobre temas de salud y hábitos de las personas en la escuela secundaria que pueden impactar en su salud.

rm( list=ls() )  #remove all objects
gc()       

          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells  541512   29    1212557 64.8         NA   669240 35.8
Vcells 1037625    8    8388608 64.0      65536  4787859 36.6

options(scipen=999)

# Carga de librerías
library(corrr)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(tidymodels)
library(rsample)
library(gridExtra)
library(knitr)
library(kableExtra)
library(GGally)

Comenzamos leyendo los datos y viendo su estructura, para esto utilizamos la función glimpse

#Importamos los datasets con los que trabajaremos
encuesta_salud_train = read.csv("/Users/dfontenla/Maestria/2022C2/EEA/practica/repo/EEA-2022/TP1/encuesta_salud_train.csv", encoding = "UTF-8") %>% mutate(id = 1:nrow(.)) 
encuesta_salud_test = read.csv("/Users/dfontenla/Maestria/2022C2/EEA/practica/repo/EEA-2022/TP1/encuesta_salud_test.csv", encoding = "UTF-8") %>% mutate(id = 1:nrow(.)) 

#encuesta_salud_train$nivel_educativo

#Observamos la estructura y variables
encuesta_salud_train %>% glimpse()

Rows: 7,024
Columns: 17
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últ…
$ id                            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…

Existen 7024 registros, compuesto de 17 datos diferentes con los que trabajaremos para entrenar diferentes modelos y poder realizar la prediccion de la variable peso. Nos quedamos con todo el dataset, como se menciona en el enunciado, corresponden a un recorte (muestra) del dataset original, luego del tratamiento de valores atípicos e ingeniería de atributos.

Análisis exploratorios

Valores únicos y porcentaje de faltantes

tabla_exploratorios =  encuesta_salud_train %>%
                                      gather(., 
                                            key = "variables", 
                                            value = "valores") %>% # agrupamos por las variables del set
                                      group_by(variables) %>% 
                                      summarise(valores_unicos = n_distinct(valores),
                                      porcentaje_faltantes = sum(is.na(valores))/nrow(encuesta_salud_train)*100) %>% 
                                      arrange(desc(porcentaje_faltantes), valores_unicos) # ordenamos por porcentaje de faltantes y valores unicos
tabla_exploratorios

Observamos un dataset que ya esta curado de datos faltantes y su composición esta dada por una mayor cantidad de datos categóricos que continuos

Análisis descriptivo

Realizamos un primer análisis descriptivo de nuestro dataset de columnas relevantes. Para ello utilizamos la función ggpairs sobre las variables numéricas con una apertura por la variable de sexo.

#Observamos la correlacion de variables
#encuesta_salud_train_numeric = encuesta_salud_train %>% dplyr::select(where(is.numeric))

# Seleccionamos las variables de interés
encuesta_salud_train_gender_numeric = encuesta_salud_train %>%
                              dplyr::select(edad, genero, altura, peso,dias_consumo_comida_rapida, consumo_diario_alcohol, 
                              dias_actividad_fisica_semanal)

encuesta_salud_train_gender_numeric %>% glimpse()

Rows: 7,024
Columns: 7
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…

encuesta_salud_train_gender_numeric %>% ggpairs(aes(color=genero), upper = list(continuous = wrap("cor", size = 3, hjust=0.8, alignPercent=0.15)), legend = 25) + 
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")

Realizamos el análisis solo numérico de la correlación entre variables

encuesta_salud_train %>% 
 correlate() %>% # convierte la matriz de corr en dataframe
  shave() %>% # solo muestra información debajo de la diagonal principal
  fashion() # acomoda los datos en forma tidy (por ej. redondeo de decimales)

Lo graficamos para tener las relaciones de una forma mas descriptiva

encuesta_salud_train %>% 
 correlate() %>% 
  network_plot(min_cor = 0.1)

Respecto a nuestra variable a predecir, el peso, observamos:

• La variable altura es aquella que tiene mayor correlación con nuestra variable a predecir
• Gráficamente también se puede observar la correlación entre altura y peso
• Exceptuando la altura, el resto de las variables presentan una correlación muy baja
• En una primera observación no parecieramos tener en el dataset la presencia de outliers que nos pueda traer complicaciones en análisis posteriores
• A pesar que existen diferencias, no observamos distribuciones significativamente diferentes según el género

Análisis frecuencia de hambre mensual

Para las categorías de la variable frecuencia de hambre mensual, analice gráficamente la distribución en términos de frecuencia relativa de:

• 1. El consumo semanal de verdura.
• 2. El consumo semanal de comida grasa.

Agrupamos los valores para obtener la frecuencia tanto de consumo semanal de verduras como de consumo semanal de comida grasa

#encuesta_salud_train_numeric <- data.frame(encuesta_salud_train)       

consumo_verduras_freq = encuesta_salud_train %>%
  group_by(frecuencia_hambre_mensual, consumo_semanal_verdura) %>%
  summarise(n = n()) %>%
  mutate(freq = n / sum(n))

consumo_verduras_freq %>% glimpse

Rows: 45
Columns: 4
Groups: frecuencia_hambre_mensual [6]
$ frecuencia_hambre_mensual <chr> "Algunas veces", "Algunas veces", "Algunas v…
$ consumo_semanal_verdura   <chr> "1 a 3 veces durante los últimos 7 días", "1…
$ n                         <int> 175, 92, 79, 21, 119, 27, 3, 73, 23, 14, 5, …
$ freq                      <dbl> 0.297113752, 0.156196944, 0.134125637, 0.035…

consumo_semanal_comida_grasa_freq = encuesta_salud_train %>%
  group_by(frecuencia_hambre_mensual, consumo_semanal_comida_grasa) %>%
  summarise(n = n()) %>%
  mutate(freq = n / sum(n))

consumo_semanal_comida_grasa_freq %>% glimpse

Rows: 43
Columns: 4
Groups: frecuencia_hambre_mensual [6]
$ frecuencia_hambre_mensual    <chr> "Algunas veces", "Algunas veces", "Alguna…
$ consumo_semanal_comida_grasa <chr> "1 a 3 veces durante los últimos 7 días",…
$ n                            <int> 252, 62, 33, 18, 81, 24, 4, 115, 27, 6, 9…
$ freq                         <dbl> 0.427843803, 0.105263158, 0.056027165, 0.…

Frecuencia de hambre mensual / consumo semanal de verduras

ggplot(data = consumo_verduras_freq, aes(x = frecuencia_hambre_mensual, y = freq, fill = consumo_semanal_verdura)) + 
geom_bar(stat = "identity", position = position_dodge()) + 
xlab("Frecuencia de hambre mensual") + 
ylab("Consumo semanal de verduras") + 
theme(axis.text.x = element_text(angle = 90)) + 
coord_flip() + 
theme_bw() + 
labs(title = "Frecuencia de hambre mensual") +
scale_fill_brewer(palette = "Paired") + 
theme(axis.text.x = element_text(size = 10)) + 
theme(axis.title.x = element_text(size = 15)) + 
theme(axis.title.y = element_text(size = 15)) + 
theme(plot.title = element_text(size = 20))

Se puede observar como para las categorias Siempre y Casi siempre de hambre mensual, aumenta de forma significativa la cantidad de valores en consumo semanal verduras, que hacen referencia a no ingerir estos alimentos con recurrencia, y si aquellas que referencian a valores semanales (1 a 3 veces durante los ultimos 7 dias y 4 a 6 veces durante los ultimos siete dias) en lugar de valores diarios

Frecuencia de hambre mensual / consumo semanal de comida grasa

ggplot(data = consumo_semanal_comida_grasa_freq, aes(x = frecuencia_hambre_mensual, y = freq, fill = consumo_semanal_comida_grasa)) + 
geom_bar(stat = "identity", position = position_dodge()) + 
xlab("Frecuencia de hambre mensual") + 
ylab("Consumo de comida grasa") + 
theme(axis.text.x = element_text(angle = 90)) + 
coord_flip() + 
theme_bw() + 
labs(title = "Frecuencia de hambre mensual") +
scale_fill_brewer(palette = "Paired") + 
theme(axis.text.x = element_text(size = 10)) + 
theme(axis.title.x = element_text(size = 15)) + 
theme(axis.title.y = element_text(size = 15)) + 
theme(plot.title = element_text(size = 20))

Por otro lado, en cuanto al análisis de la relación entre el consumo de comida grasa y la frecuencia de hambre, observamos que el consumo habitual de comidad grasa, aumenta la frecuencia de hambre mensual. Buscando referencias a esta conclusión, existen estudios relacionados a la obecidad donde focalizan sobre esta relación donde los alimentos con altos contenidos grasos, dan más hambre

Modelo inicial

Creación del modelo

Se plantea que una primera alternativa para modelar el peso es: E(peso) = β0 +β1altura+β2edad+β3genero+β4diasActividadF isicaSemanal+β5consumoDiarioAlcohol

Realizamos una selección de las variables de nuestro interés para poder realizar el modelado, y observamos su estructura

#encuesta_salud_train %>% glimpse
encuesta_salud_modelo_base = encuesta_salud_train %>%
                              # Seleccionamos las variables de interés
                              dplyr::select(edad, genero, altura, peso,dias_actividad_fisica_semanal,consumo_diario_alcohol)

encuesta_salud_modelo_base %>% glimpse

Rows: 7,024
Columns: 6
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …

Vemos cómo es la correlación entre las variables numéricas.

encuesta_salud_modelo_base %>% ggpairs(aes(color=genero), upper = list(continuous = wrap("cor", size = 3, hjust=0.8, alignPercent=0.15)), legend = 25) + 
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")

Armamos un modelo para predecir el peso en función de la altura, edad, género, días de actividad física semanal y consumo diario de alcohol.

E(peso) = β0 +β1altura+β2edad+β3genero+β4diasActividadFisicaSemanal+β5consumoDiarioAlcohol

# ajustamos modelo lineal múltiple
modelo_ln_base <- lm(peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol, data = encuesta_salud_modelo_base)
# Resumen del modelo
tidy_ln_base <- tidy(modelo_ln_base, conf.int = TRUE)
tidy_ln_base

Evaluación de coeficientes

Significado de los coeficientes estimados

• El valor de β0^ (ordenada al origen) es -68.92 kg, lo que corresponde al peso esperado de una persona femenina sin edad, sin altura, sin dias de actividad fisica y sin consumo diario en alcohol. Lo cual, este caso carecería de sentido ya que las personas femeninas deberían tener como valores de altura y edad numeros positivos > 0.

• β0^ + βgeneroMasculino es la media del peso para las personas de género masculino, dada la edad, altura, dias de consumo de alcohol, dias de actividad fisica. Por lo tanto, βgeneroMasculino es la diferencia en los niveles medios de pesos de las personas masculinas respecto de las femeninas (categoría basal). Es decir, βgeneroMasculino (1.262643558) indica cuánto más alta es la función de respuesta (peso) para las personas masculinas respecto de las femeninas (categoría basal), dada la edad, altura, dias de consumo de alcohol y dias de actividad física.

• El coeficiente estimado de edad es de 1.4kg, lo que indica que si mantenemos el número de altura, género, dias de actividad fisica y consumo diario de alchol, cada incremento adicional de edad corresponde a un aumento de 1.4kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma altura, dias de actividad fisica, genero y consumo de alcohol, pero teniendo una un año más de edad que la otra, el peso esperado para la de mayor edad será 1.4kg más que la de menor pesaje.

• El coeficiente estimado de altura es de 0.65kg, lo que indica que si mantenemos el número de género, edad, dias de actividad física y consumo diario de alchol, cada incremento adicional de altura corresponde a un aumento de 0.65kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma edad, dias de actividad fisica, genero y consumo de alcohol, pero teniendo una un centímetro más de altura que la otra, el peso esperado para la de mayor altura será 0.65kg más que la de menor pesaje.

• El coeficiente estimado de dias de actividad fisica es de -0.0874kg, lo que indica que si mantenemos el género, edad, altura y consumo diario de alchol, cada incremento de un día de actividad fisica adicional corresponde a un decremento de 0.0874kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma edad, altura, género y consumo de alcohol, pero teniendo una un día de actividad fisica más que la otra, el peso esperado para la de mayor dias de actividad física será 0.0874kg menor que la de mayor pesaje.

• El coeficiente estimado de dias de consumo de alcohol es de 0.00727kg, lo que indica que si mantenemos el número de género, edad, altura y días de actividad física, cada incremento de un día de consumo de alcohol adicional corresponde a un aumento de 0.00727kg, en promedio en el peso de la persona femenina. O lo que es igual, dadas dos personas con la misma edad, altura, género y días de actividad física, pero teniendo una un día de consumo de alcohol más que la otra, el peso esperado para la de mayor días de consumo de alcohol será 0.00727kg más que la de menor pesaje.

Significancia individual y global

Inferencia de los βk (test de significatividad individual)

Para evaluar la significatividad individual de cada una de las variables se analiza el test t que busca probar si el coeficiente de regresión correspondiente a dicha variable es distinto de 0 (figura en la tabla resumen de resultados de la regresión).

Es decir, buscamos probar:

H0:βk=0 H1:βk≠0.

options("scipen"=1)
tidy_ln_base %>% select(term, statistic, p.value, conf.low, conf.high)

significancia = tidy_ln_base %>% dplyr::select(term, statistic, p.value, conf.low, conf.high)

significancia

Evaluando la significancia de cada variable observamos lo siguiente

Para las variables edad, altura, generoMasculino, generoFemenino (basal):

• Las variables resultan estadísticamente significativas para explicar el peso de las personas (p-valores < 0.05).
• Además del resultado del test, podemos corroborar que los intervalos de confianza del 95% para los coeficientes estimados no contienen al 0 en ninguno de los casos.

Para las variables días de actividad física semanal y consumo diario de alcohol:

• Las variables NO resultan estadísticamente significativas para explicar el peso de las personas (p-valores > 0.05).
• Además del resultado del test, podemos corroborar que los intervalos de confianza del 95% para los coeficientes estimados contienen al 0 en los casos.

tidy(anova(modelo_ln_base))

La tabla de ANOVA muestra que, según el resultado del test F, la variable género en su conjunto, también resulta estadísticamente significativa para explicar al precio (p-valor < 0.05).

Test F (test de significatividad global)

Se construye para testear las hipótesis:

H0:β1=β2=···=βp−1=0

H1: no todos los βk (k=1,2,…,p−1) son iguales a 0.

Observemos que H0 dice que no hay vínculo entre la variable respuesta y las regresoras. En cambio, H1 dice que al menos una de las variables regresoras sirve para predecir a Y. Los resultados de este test se pueden observar haciendo un summary() del modelo o glance().

glance(modelo_ln_base)

F-statistic: 770.4 on 5 and 7018 DF, p-value: < 2.2e-16

Test significativo por lo que determinamos que se rechaza la hipotesis nula, nuestro modelo es relevante en comparación al modelo base, para determinar que arroja resultados significativos tomando la variable predictora

El R-cuadrado permite medir el porcentaje de variabilidad del fenómeno que el modelo logra explicar. Por este motivo es una métrica que nos permite evaluar la capacidad explicativa del modelo y poder comparar modelos entre sí bajo ciertas condiciones.

En nuestro caso tenemos un valor de Multiple R-squared: 0.3544, Adjusted R-squared: 0.3539

Modelo categóricos

Creación modelo categórico

Se sugiere probar un modelo que incopore el consumo semanal de snacks y una interacción entre el género y la edad, en lugar de actividad física y consumo de alcohol:

E(peso) = β0 + β1altura + β2edad + β3genero + β4consumoSemanalSnacks + β5genero · edad

encuesta_salud_train %>% glimpse

Rows: 7,024
Columns: 17
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últ…
$ id                            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…

encuesta_salud_second_model = encuesta_salud_train %>%
                              # Seleccionamos las variables de interés
                              dplyr::select(edad, genero, altura, peso,consumo_semanal_snacks)

encuesta_salud_second_model %>% glimpse

Rows: 7,024
Columns: 5
$ edad                   <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16, 14,…
$ genero                 <chr> "Femenino", "Masculino", "Masculino", "Masculin…
$ altura                 <int> 165, 178, 172, 170, 170, 178, 156, 163, 164, 16…
$ peso                   <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100, 33…
$ consumo_semanal_snacks <chr> "1 a 3 veces durante los últimos 7 días", "No c…

Dejamos al valor No comí comida salada o snacks en los últimos 7 días de consumo semanal de snacks como categoría basal en una nueva variable que llamaremos consumo_semanal_snacks_relevel

encuesta_salud_second_model$consumo_semanal_snacks_relevel <- relevel(as.factor(encuesta_salud_second_model$consumo_semanal_snacks), ref = 8)

Graficámos la variable categórica consumoSemanalSnacks para analizar sus valores en función del peso

# armo boxplots paralelos consumo de snacks en función del peso
ggplot( encuesta_salud_second_model, aes(x = fct_reorder(consumo_semanal_snacks_relevel, peso, .desc = T), y = peso)) + 
  geom_boxplot(alpha = 0.75, aes(fill = consumo_semanal_snacks_relevel)) + 
  theme_minimal() + 
  theme(legend.position = 'none')+
  labs(y = "Peso", x = "Consumo snacks")  +
  ggtitle("Boxplots del peso en función del consumo de snacks")+
  theme (axis.text.x = element_text(face="italic", colour="dark grey", size = 8, angle = 90))

modelo_ln_2_r <- lm(peso ~ consumo_semanal_snacks_relevel + edad + genero + altura + (edad * genero), data = encuesta_salud_second_model)

# Resumen del modelo
tidy_ln_2_r <- tidy(modelo_ln_2_r, conf.int = TRUE)
print(tidy_ln_2_r)

# A tibble: 12 × 7
   term                estimate std.error statistic   p.value conf.low conf.high
   <chr>                  <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
 1 (Intercept)          -64.2      2.83      -22.7  4.08e-110 -69.7      -58.7  
 2 consumo_semanal_sn…   -1.35     0.276      -4.89 1.03e-  6  -1.89      -0.810
 3 consumo_semanal_sn…   -0.608    0.456      -1.33 1.83e-  1  -1.50       0.286
 4 consumo_semanal_sn…   -1.09     0.685      -1.60 1.10e-  1  -2.44       0.247
 5 consumo_semanal_sn…   -1.28     1.01       -1.26 2.07e-  1  -3.26       0.707
 6 consumo_semanal_sn…   -2.27     0.450      -5.05 4.64e-  7  -3.15      -1.39 
 7 consumo_semanal_sn…   -2.57     0.881      -2.91 3.59e-  3  -4.29      -0.840
 8 consumo_semanal_sn…   -4.44     1.95       -2.28 2.29e-  2  -8.27      -0.615
 9 edad                   1.22     0.121      10.1  7.55e- 24   0.986      1.46 
10 generoMasculino       -4.61     2.68       -1.72 8.52e-  2  -9.86       0.640
11 altura                 0.643    0.0145     44.3  0           0.614      0.671
12 edad:generoMasculi…    0.391    0.179       2.19 2.88e-  2   0.0405     0.742

Interpretación de coeficientes / Análisis de significancia

Intercept

• El valor de β0 (ordenada al origen) es -64.2 kg, lo que corresponde al peso esperado de una persona femenina, que no comió comida salada o snacks en los últimos 7 dias y contiene la suma del coeficiente edad:generoFemenino manteniendo el resto de las variables constantes. (Una variable de cada atributo categoróco forman parte del Intercept, en este caso género femenino, aquella que redefinimos para que sea basal de consumo de snacks No comió comida salada o snacks en los últimos 7 dias y la nueva variable creada edad:generoFemenino)

Coeficiente edad:generoMasculino:

• El coeficiente estimado de edad:generoMasculino es de 0.391kg, lo que indica que si mantenemos el número de altura, genero, edad y consumo semanal de snacks, cada incremento de una unidad para edad:generoMasculino corresponde a un aumento de 0.391kg, en promedio en el peso de la persona femenina. => Significativa

Coeficiente del consumo semanal de snacks:

• El coeficiente estimado de consumo semanal snacks 1 a 3 veces durante los últimos 7 días es de -1.35, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 1 a 3 veces durante los últimos 7 días, corresponde a una baja de 1.35kg, en promedio en el peso de la persona. => Significativa

• El coeficiente estimado de consumo_semanal_snacks 1 vez al día es de -0.608, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 1 vez al día corresponde a una baja de 0.608, en promedio en el peso de la persona. No significativa

• El coeficiente estimado de consumo_semanal_snacks 2 veces al día es de -1.09, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 2 veces al día corresponde a una baja de 1.09, en promedio en el peso de la persona. => No significativa

• El coeficiente estimado de consumo_semanal_snacks 3 veces al día es de -1.28, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 3 veces al día corresponde a una baja de 1.28, en promedio en el peso de la persona. => No significativa

• El coeficiente estimado de consumo_semanal_snacks 4 a 6 veces durante los últimos 7 días es de -2.27, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacks 4 a 6 veces durante los últimos 7 días corresponde a una baja de 2.27kg, en promedio en el peso de la persona. => Significativa

• El coeficiente estimado de consumo_semanal_snacks 4 o más veces al día es de -2.57, si mantenemos el resto de las variables constantes, el incremento de una unidad para estimado de consumo_semanal_snacks 4 o más veces al día corresponde a una baja de 2.57kg, en promedio en el peso de la persona. => Significativa

• El coeficiente estimado de consumo_semanal_snacks Dato perdido es de -4.44, si mantenemos el resto de las variables constantes, el incremento de una unidad para consumo_semanal_snacksDato perdido corresponde a una baja de 4.44kg, en promedio en el peso de la persona. => Significativa

Los coeficientes para las diferentes categorias de consumo semanal de snacks algunos son significativos y otros no.

Realizamos el test de anova para ver la significancia de las variables en su conjunto

tidy(anova(modelo_ln_2_r))

La tabla de ANOVA muestra que, según el resultado del test F, la variable consumo_semanal_snacks en su conjunto resulta estadísticamente significativa para explicar al precio (p-valor < 0.05).

Redefinicion de categóricas

Creamos una variable nueva que agrupe las tres variables categóricas que nos dieron una NO significancia en el análisis previo, las mismas serian

• [3] “consumo_semanal_snacks_relevel1 vez al día”

• [4] “consumo_semanal_snacks_relevel2 veces al día”

• [5] “consumo_semanal_snacks_relevel3 veces al día”

encuesta_salud_train %>% glimpse

Rows: 7,024
Columns: 17
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últ…
$ id                            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…

encuesta_salud_second_model_agrupado = encuesta_salud_train %>%
                              # Seleccionamos las variables de interés
                              dplyr::select(edad, genero, altura, peso,consumo_semanal_snacks)

encuesta_salud_second_model_agrupado$consumo_semanal_snacks_relevel <- relevel(as.factor(encuesta_salud_second_model_agrupado$consumo_semanal_snacks), ref = 8)

Redefinimos la variable agrupando los valores categorizados como 1 vez al día, 2 veces al día y tres veces al día, en 1 a 3 veces al dia.

encuesta_salud_second_model_agrupado = encuesta_salud_second_model_agrupado %>% mutate(consumo_semanal_snacks = case_when(startsWith(consumo_semanal_snacks, "1 vez al día") ~ "1 a 3 veces al dia",
                                 startsWith(consumo_semanal_snacks, "2 veces al día") ~ "consumo_semanal_snacks 1 a 3 veces al dia",
                                 startsWith(consumo_semanal_snacks, "3 veces al día") ~ "consumo_semanal_snacks 1 a 3 veces al dia",
                                 TRUE ~ consumo_semanal_snacks))

#encuesta_salud_second_model_agrupado$consumo_semanal_snacks

modelo_ln_2_r_agrupado <- lm(peso ~ consumo_semanal_snacks + edad + genero + altura + (edad * genero), data = encuesta_salud_second_model_agrupado)

# Resumen del modelo
tidy_ln_2_r_agrupado <- tidy(modelo_ln_2_r_agrupado, conf.int = TRUE)
print(tidy_ln_2_r_agrupado)

# A tibble: 11 × 7
   term                estimate std.error statistic   p.value conf.low conf.high
   <chr>                  <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
 1 (Intercept)          -64.8      2.83     -22.9   4.14e-112 -70.4     -59.3   
 2 consumo_semanal_sn…   -0.744    0.440     -1.69  9.11e-  2  -1.61      0.119 
 3 consumo_semanal_sn…   -1.66     0.565     -2.94  3.28e-  3  -2.77     -0.554 
 4 consumo_semanal_sn…   -1.96     0.944     -2.07  3.81e-  2  -3.81     -0.108 
 5 consumo_semanal_sn…   -0.542    0.676     -0.801 4.23e-  1  -1.87      0.784 
 6 consumo_semanal_sn…   -3.83     1.98      -1.93  5.30e-  2  -7.72      0.0503
 7 consumo_semanal_sn…    0.608    0.456      1.33  1.83e-  1  -0.286     1.50  
 8 edad                   1.22     0.121     10.1   7.57e- 24   0.986     1.46  
 9 generoMasculino       -4.61     2.68      -1.72  8.52e-  2  -9.86      0.640 
10 altura                 0.643    0.0145    44.3   0           0.614     0.671 
11 edad:generoMasculi…    0.391    0.179      2.19  2.88e-  2   0.0405    0.742 

Observamos que la redefinición de las variables no significativas generando un agrupamiento para de las mismas mantiene la no significancia de estos valores para predecir el peso de una persona.

Evaluación de modelo

Ejecutamos las funcines tidy, glance y summary para obtener tanto la significacia de los atributos del modelo, la significancia del modelo en si (Test F) y la variabilidad del mismo (R-cuadrado y R-cuadrado ajustado).

options("scipen"=1)
tidy_ln_2_r %>%
  dplyr::select(term, statistic, p.value, conf.low, conf.high)

glance(modelo_ln_2_r)

summary(modelo_ln_2_r)

Call:
lm(formula = peso ~ consumo_semanal_snacks_relevel + edad + genero + 
    altura + (edad * genero), data = encuesta_salud_second_model)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.005  -6.452  -1.157   5.052  75.767 

Coefficients:
                                                                      Estimate
(Intercept)                                                          -64.19940
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días  -1.35163
consumo_semanal_snacks_relevel1 vez al día                            -0.60788
consumo_semanal_snacks_relevel2 veces al día                          -1.09466
consumo_semanal_snacks_relevel3 veces al día                          -1.27596
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días  -2.27004
consumo_semanal_snacks_relevel4 o más veces al día                    -2.56697
consumo_semanal_snacks_relevelDato perdido                            -4.44042
edad                                                                   1.22309
generoMasculino                                                       -4.61150
altura                                                                 0.64292
edad:generoMasculino                                                   0.39146
                                                                     Std. Error
(Intercept)                                                             2.82849
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días    0.27645
consumo_semanal_snacks_relevel1 vez al día                              0.45601
consumo_semanal_snacks_relevel2 veces al día                            0.68457
consumo_semanal_snacks_relevel3 veces al día                            1.01150
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días    0.44992
consumo_semanal_snacks_relevel4 o más veces al día                      0.88112
consumo_semanal_snacks_relevelDato perdido                              1.95125
edad                                                                    0.12102
generoMasculino                                                         2.67877
altura                                                                  0.01453
edad:generoMasculino                                                    0.17902
                                                                     t value
(Intercept)                                                          -22.697
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días  -4.889
consumo_semanal_snacks_relevel1 vez al día                            -1.333
consumo_semanal_snacks_relevel2 veces al día                          -1.599
consumo_semanal_snacks_relevel3 veces al día                          -1.261
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días  -5.045
consumo_semanal_snacks_relevel4 o más veces al día                    -2.913
consumo_semanal_snacks_relevelDato perdido                            -2.276
edad                                                                  10.106
generoMasculino                                                       -1.722
altura                                                                44.250
edad:generoMasculino                                                   2.187
                                                                     Pr(>|t|)
(Intercept)                                                           < 2e-16
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días 1.03e-06
consumo_semanal_snacks_relevel1 vez al día                            0.18257
consumo_semanal_snacks_relevel2 veces al día                          0.10985
consumo_semanal_snacks_relevel3 veces al día                          0.20719
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días 4.64e-07
consumo_semanal_snacks_relevel4 o más veces al día                    0.00359
consumo_semanal_snacks_relevelDato perdido                            0.02289
edad                                                                  < 2e-16
generoMasculino                                                       0.08520
altura                                                                < 2e-16
edad:generoMasculino                                                  0.02880
                                                                        
(Intercept)                                                          ***
consumo_semanal_snacks_relevel1 a 3 veces durante los últimos 7 días ***
consumo_semanal_snacks_relevel1 vez al día                              
consumo_semanal_snacks_relevel2 veces al día                            
consumo_semanal_snacks_relevel3 veces al día                            
consumo_semanal_snacks_relevel4 a 6 veces durante los últimos 7 días ***
consumo_semanal_snacks_relevel4 o más veces al día                   ** 
consumo_semanal_snacks_relevelDato perdido                           *  
edad                                                                 ***
generoMasculino                                                      .  
altura                                                               ***
edad:generoMasculino                                                 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.887 on 7012 degrees of freedom
Multiple R-squared:  0.3585,    Adjusted R-squared:  0.3575 
F-statistic: 356.3 on 11 and 7012 DF,  p-value: < 2.2e-16

El porcentaje de variabilidad que explica este nuevo modelo es de R-cuadrado = 0.359, R-cuadrado ajustado = 0.358 Evaluando la significativida del modelo, obtenemos F-statistic: 356.3 on 11 and 7012 DF, p-value: < 2.2e-16 Modelo significativo

Modelos propios y evaluación

Creación de nuevos modelos

Proponemos un modelo que evalue el peso a partir del nivel educativo, consumo_semanal_comida_grasa, edad, genero, altura y una interrelación de altura, edad y género.

f(peso) ~ nivel_educativo + consumo_semanal_comida_grasa + edad + genero + altura + (altura * edad * genero)

Seleccionamos nuestras variables de interés

encuesta_salud_train %>% glimpse

Rows: 7,024
Columns: 17
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últ…
$ id                            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…

encuesta_salud_third_model = encuesta_salud_train %>%
                              # Seleccionamos las variables de interés
                              dplyr::select(edad, genero, altura, peso,nivel_educativo,consumo_semanal_comida_grasa)

encuesta_salud_third_model %>% glimpse

Rows: 7,024
Columns: 6
$ edad                         <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 1…
$ genero                       <chr> "Femenino", "Masculino", "Masculino", "Ma…
$ altura                       <int> 165, 178, 172, 170, 170, 178, 156, 163, 1…
$ peso                         <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 1…
$ nivel_educativo              <chr> "2do año/11vo grado nivel Polimodal o 4to…
$ consumo_semanal_comida_grasa <chr> "No comí comida alta en grasa en los últi…

Generamos el modelo lineal

modelo_ln_3_r <- lm(peso ~ nivel_educativo + consumo_semanal_comida_grasa + edad + genero + altura + (altura * edad * genero), data = encuesta_salud_third_model)
# Resumen del modelo
tidy_ln_3_r <- tidy(modelo_ln_3_r, conf.int = TRUE)
print(tidy_ln_3_r)

# A tibble: 20 × 7
   term                  estimate std.error statistic p.value conf.low conf.high
   <chr>                    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
 1 (Intercept)           -32.9      38.6      -0.852  0.394   -1.09e+2   42.8   
 2 nivel_educativo2do a…  -0.779     0.356    -2.19   0.0287  -1.48e+0   -0.0808
 3 nivel_educativo3er a…  -0.423     0.411    -1.03   0.304   -1.23e+0    0.383 
 4 nivel_educativo8vo g…  -0.930     0.458    -2.03   0.0424  -1.83e+0   -0.0321
 5 nivel_educativo9no g…  -0.400     0.379    -1.06   0.290   -1.14e+0    0.342 
 6 nivel_educativoDato …  -1.38      1.02     -1.35   0.176   -3.39e+0    0.622 
 7 consumo_semanal_comi…   1.10      0.428     2.57   0.0102   2.61e-1    1.94  
 8 consumo_semanal_comi…  -1.42      0.630    -2.26   0.0239  -2.66e+0   -0.188 
 9 consumo_semanal_comi…  -1.33      0.986    -1.35   0.177   -3.27e+0    0.601 
10 consumo_semanal_comi…  -0.915     0.349    -2.62   0.00870 -1.60e+0   -0.231 
11 consumo_semanal_comi…  -0.0676    0.849    -0.0796 0.937   -1.73e+0    1.60  
12 consumo_semanal_comi…  -2.94      1.40     -2.10   0.0354  -5.68e+0   -0.200 
13 consumo_semanal_comi…   0.0307    0.308     0.0996 0.921   -5.74e-1    0.635 
14 edad                   -0.343     2.57     -0.133  0.894   -5.38e+0    4.70  
15 generoMasculino        11.1      50.3       0.221  0.825   -8.75e+1  110.    
16 altura                  0.440     0.241     1.82   0.0682  -3.29e-2    0.913 
17 edad:altura             0.0101    0.0161    0.629  0.529   -2.14e-2    0.0416
18 generoMasculino:altu…  -0.0673    0.308    -0.219  0.827   -6.71e-1    0.536 
19 edad:generoMasculino   -1.74      3.38     -0.516  0.606   -8.37e+0    4.88  
20 edad:generoMasculino…   0.0111    0.0206    0.538  0.591   -2.93e-2    0.0515

glance(modelo_ln_3_r)

tidy(anova(modelo_ln_3_r))

Proponemos un modelo que evalue el peso a partir del nivel educativo, consumo_semanal_comida_grasa, genero, altura y una interrelación de dias consumo comida rápida y genero f(peso) ~ nivel_educativo + consumo_semanal_comida_grasa + genero + altura + (dias_consumo_comida_rapida * genero)

Seleccionamos nuestras variables de interés

encuesta_salud_train %>% glimpse

Rows: 7,024
Columns: 17
$ record                        <int> 502, 26488, 31473, 14154, 36578, 53730, …
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo               <chr> "2do año/11vo grado nivel Polimodal o 4t…
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, …
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, …
$ frecuencia_hambre_mensual     <chr> "Rara vez", "Rara vez", "Nunca", "Nunca"…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1…
$ edad_consumo_alcohol          <chr> "14 o 15 años", "7 años o menos", "Nunca…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, …
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1…
$ consumo_semanal_frutas        <chr> "No comí frutas durante los últimos 7 dí…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas      <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa  <chr> "No comí comida alta en grasa en los últ…
$ id                            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…

encuesta_salud_fourth_model = encuesta_salud_train %>%
                              # Seleccionamos las variables de interés
                              dplyr::select(edad, genero, altura, peso,nivel_educativo,dias_consumo_comida_rapida)

encuesta_salud_fourth_model %>% glimpse

Rows: 7,024
Columns: 6
$ edad                       <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16,…
$ genero                     <chr> "Femenino", "Masculino", "Masculino", "Masc…
$ altura                     <int> 165, 178, 172, 170, 170, 178, 156, 163, 164…
$ peso                       <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100…
$ nivel_educativo            <chr> "2do año/11vo grado nivel Polimodal o 4to a…
$ dias_consumo_comida_rapida <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1, 3…

Generamos nuestro modelo lineal

modelo_ln_4_r <- lm(peso ~ nivel_educativo + dias_consumo_comida_rapida + genero + altura + (dias_consumo_comida_rapida * genero), data = encuesta_salud_fourth_model)
# Resumen del modelo
tidy_ln_4_r <- tidy(modelo_ln_4_r, conf.int = TRUE)
print(tidy_ln_4_r)

# A tibble: 10 × 7
   term                estimate std.error statistic   p.value conf.low conf.high
   <chr>                  <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
 1 (Intercept)         -5.12e+1    2.31    -22.2    3.15e-105  -55.7    -46.6   
 2 nivel_educativo2do…  1.86e-1    0.348     0.534  5.93e-  1   -0.496    0.867 
 3 nivel_educativo3er…  1.50e+0    0.372     4.03   5.52e-  5    0.771    2.23  
 4 nivel_educativo8vo… -2.76e+0    0.425    -6.48   1.00e- 10   -3.59    -1.92  
 5 nivel_educativo9no… -1.52e+0    0.370    -4.10   4.20e-  5   -2.24    -0.791 
 6 nivel_educativoDat… -8.26e-1    1.03     -0.801  4.23e-  1   -2.85     1.20  
 7 dias_consumo_comid… -5.82e-3    0.117    -0.0498 9.60e-  1   -0.235    0.223 
 8 generoMasculino      1.41e+0    0.306     4.62   3.86e-  6    0.815    2.01  
 9 altura               6.72e-1    0.0143   47.0    0            0.644    0.700 
10 dias_consumo_comid… -4.14e-1    0.175    -2.37   1.79e-  2   -0.757   -0.0712

glance(modelo_ln_4_r)

tidy(anova(modelo_ln_4_r))

Evaluación de modelos

Agrupamos todos los modelos que nos interesa comparar, en este caso es el modelo inicial, el modelo categóricas con las categorías redefinidas de la variable consumoSemanalSnacks y los dos modelos desarrollados en este punto, interrelación (altura * edad * genero) y (dias_consumo_comida_rapida * genero)

# armamos lista con todos los modelos
models <- list(modelo_ln_4_r = modelo_ln_4_r, modelo_ln_3_r = modelo_ln_3_r, modelo_ln_2_r = modelo_ln_2_r, modelo_ln_base = modelo_ln_base)
# calculamos las variables resumen
map_df(models, tidy, .id = "model")

Selección del mejor modelo en training

Utilizamos la función augment para predecir el peso sobre el dataset de training

lista_predicciones_training = map(.x = models, .f = augment) 

Calculamos el RMSE para todos los modelos a evaluar

map_dfr(.x = lista_predicciones_training, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Calculamos el MAE para todos los modelos a evaluar

map_dfr(.x = lista_predicciones_training, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Obtenemos el R-cuadrado y R-cuadrado ajustado para poder comparar la variabilidad de los modelos

df_evaluacion_train = map_df(models, glance, .id = "model") %>%
  # ordenamos por R2 ajustado
  arrange(desc(adj.r.squared))

df_evaluacion_train

Concluimos sobre el dataset de training

• Analizando el R-cuadrado ajustado, el mejor modelo resulta ser modelo_ln_2_r que explica el peso en función de la altura, edad, género, dias de actividad física semanal y consumo diario de alcohol
• Este modelo (modelo_ln_2_r) explica 35,8% de la variabilidad del peso, es decir, más que todos los restantes modelos.
• El modelo con menor RMSE es el modelo_ln_3_r con un valor de 9.87, no coincide con el mejor R-cuadrado ajustado, pero si por una diferencia mínima modelo_ln_2_r (9.88)
• En cuanto al MAE, tanto el modelo modelo_ln_2_r como modelo_ln_3_r, son los que tienen mejores resultados con 7.43

Evaluación en el dataset de testing

Utilizamos la función augment para predecir el peso sobre el dataset de training

# Aplicamos la función augment a los 4 modelos con el set de testing
encuesta_salud_test$consumo_semanal_snacks_relevel <- relevel(as.factor(encuesta_salud_test$consumo_semanal_snacks), ref = 8)
lista_predicciones_testing = map(.x = models, .f = augment, newdata = encuesta_salud_test) 

Calculamos el RMSE para todos los modelos a evaluar con las predicciones en training

# Obtenemos el RMSE para los 4 modelos
map_dfr(.x = lista_predicciones_testing, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Calculamos el MAE para todos los modelos a evaluar con las predicciones en training

map_dfr(.x = lista_predicciones_testing, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Concluimos sobre el dataset de testing:

• Tanto el modelo modelo_ln_2_r como el modelo_ln_3_r son los que mejor RMSE tienen como en la medición sobre training con 10.2
• Tanto el modelo modelo_ln_2_r como el modelo_ln_3_r son los que mejor MAE tienen como en la medición sobre training con 7.56

Diagnóstico del modelo

Realizaremos diferentes pruebas para validar el cumplimiento o no de los supuestos del modelo lineal, estos consisten de: εi∼N(0,σ2) independientes entre sí. Los errores son inobservables, por lo tanto tendremos que trabajar con su correlato empírico: los residuos en las técnicas de diagnóstico.

Procedemos a graficar el comportamiento de los residuos del modelo

au_modelos = map_df(models, augment, .id = "model")

g1 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"), 
       aes(.fitted, .resid)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE) +
  labs(title = "Residuos vs valores predichos MODELO BASE") + 
  theme_bw()

g2 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"), 
       aes(sample = .std.resid)) +
  stat_qq() +
  geom_abline() +
  labs(title = "Normal QQ plot MODELO BASE") + 
  theme_bw()

g3 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"), 
       aes(.fitted, sqrt(abs(.std.resid)))) +
  geom_point() +
  geom_smooth(se = FALSE) + 
  theme_bw() +
  labs(title = "Scale-location plot MODELO BASE")

g4 = ggplot(au_modelos %>% filter(model == "modelo_ln_base"), 
       aes(.hat, .std.resid)) +
  geom_vline(size = 2, colour = "white", xintercept = 0) +
  geom_hline(size = 2, colour = "white", yintercept = 0) +
  geom_point() + 
  geom_smooth(se = FALSE) + 
  theme_bw() +
  labs(title = "Residual vs leverage MODELO BASE")
# grafico todos juntos
grid.arrange(g1, g2, g3, g4, nrow = 2)

Conclusiones:

• Residuos vs valores predichos: Parece no existir estructura en los datos, la varianza no parece incrementarse con los valores predichos, por lo que satisface el supuesto de homocedasticidad.
• Normal QQ plot: Tanto el extremo superior derecho en particular, como el izquierdo inferior no se ajustan a la distribución teórica, no siguen una distribución normal
• Residual vs leverage: Existen algunos puntos de alto leverage, habria que realizar un analisis mas profundo para entender la influencia de estos valores sobre el modelo, gráficamente no es simple determinarlo

Diagnóstico del modelo: El modelo creado no cumple con los supuestos del modelo lineal. Parecen existir problemas de falta de normalidad y presencia de observaciones de alto leverage.

Modelo robusto

Análisis del dataset

Leemos el nuevo archivo con la incorporación de algunas observaciones adicionales que pueden incluir valores atípicos

library(MASS)

encuesta_salud_train_original = read_csv("/Users/dfontenla/Maestria/2022C2/EEA/practica/repo/EEA-2022/TP1/encuesta_salud_modelo6.csv") %>% mutate(id = 1:nrow(.)) 
encuesta_salud_train_original %>% glimpse()

Rows: 7,129
Columns: 17
$ record                        <dbl> 48568, 39818, 11360, 50320, 26325, 44811…
$ edad                          <dbl> 17, 14, 16, 14, 15, 13, 16, 13, 15, 14, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ nivel_educativo               <chr> "3er año/12vo grado nivel Polimodal o 5t…
$ altura                        <dbl> 160, 189, 166, 178, 167, 164, 162, 170, …
$ peso                          <dbl> 51, 75, 50, 80, 50, 55, 60, 62, 48, 71, …
$ frecuencia_hambre_mensual     <chr> "Nunca", "Nunca", "Nunca", "Nunca", "Nun…
$ dias_consumo_comida_rapida    <dbl> 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0…
$ edad_consumo_alcohol          <chr> "12 o 13 años", "Nunca tomé alcohol más …
$ consumo_diario_alcohol        <dbl> 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, …
$ dias_actividad_fisica_semanal <dbl> 2, 0, 0, 2, 1, 3, 0, 2, 3, 7, 2, 2, 1, 6…
$ consumo_semanal_frutas        <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_verdura       <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_gaseosas      <chr> "4 a 6 veces durante los últimos 7 días"…
$ consumo_semanal_snacks        <chr> "1 a 3 veces durante los últimos 7 días"…
$ consumo_semanal_comida_grasa  <chr> "4 a 6 veces durante los últimos 7 días"…
$ id                            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1…

Graficamos la relación peso / altura para observar las anomalias insertadas en el nuevo dataset

encuesta_salud_train_original %>% ggplot(., aes(x = altura, y = peso)) + 
  geom_point(alpha=0.5) + #capa de los datos
  theme_bw() +
  labs(title="Modelo Lineal Múltiple: Altura x Peso", x="Altura", y="Peso") 

Observamos que existen una cantidad de puntos relevantes que se alejan de forma categórica de la distribución general, una importante cantidad de outliers

Nos quedamos con las variables de interés para reproducir el modelo base

encuesta_salud_train_original_base = encuesta_salud_train_original %>%
                              # Seleccionamos las variables de interés
                              dplyr::select(edad, genero, altura, peso,dias_actividad_fisica_semanal,consumo_diario_alcohol)

encuesta_salud_train_original_base %>% glimpse

Rows: 7,129
Columns: 6
$ edad                          <dbl> 17, 14, 16, 14, 15, 13, 16, 13, 15, 14, …
$ genero                        <chr> "Femenino", "Masculino", "Masculino", "M…
$ altura                        <dbl> 160, 189, 166, 178, 167, 164, 162, 170, …
$ peso                          <dbl> 51, 75, 50, 80, 50, 55, 60, 62, 48, 71, …
$ dias_actividad_fisica_semanal <dbl> 2, 0, 0, 2, 1, 3, 0, 2, 3, 7, 2, 2, 1, 6…
$ consumo_diario_alcohol        <dbl> 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, …

Graficamos las relaciones entre variables del nuevo dataset

encuesta_salud_train_original_base %>% ggpairs(aes(color=genero), upper = list(continuous = wrap("cor", size = 3, hjust=0.8, alignPercent=0.15)), legend = 25) + 
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5), legend.position = "bottom")

Representamos la correlación de Pearson entre variables

encuesta_salud_train_original_base %>% 
 correlate() %>% # convierte la matriz de corr en dataframe
  shave() %>% # solo muestra información debajo de la diagonal principal
  fashion() # acomoda los datos en forma tidy (por ej. redondeo de decimales)

Podemos observar que la correlación entre peso y altura en relación al dataset base, disminuyó considerablemente pasando de un 0.5 a un 0.29, gráficamente se visualiza una clara correspondencia entre estos valores, ahora bien para el nuevo dataset existe una cantidad de datos outliers significativos

Creación del modelo

Ajustamos el dataset a nuestro modelo lineal múltiple base (peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol)

# ajustamos modelo lineal múltiple
modelo_ds_original_ln_base <- lm(peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol, data = encuesta_salud_train_original_base)
# Resumen del modelo
tidy_ln_ds_original_base <- tidy(modelo_ds_original_ln_base, conf.int = TRUE)
tidy_ln_ds_original_base

Realizamos un diagnóstico del modelo a través de los residuos

plot(modelo_ds_original_ln_base)

Vemos la significancia de los atributos

options("scipen"=1)
tidy_ln_ds_original_base %>%
  dplyr::select(term, statistic, p.value, conf.low, conf.high)

Analizamos la variabilidad explicativa del modelo

glance(modelo_ds_original_ln_base)

Observamos que la variabilidad porcentual explicada por R-cuadrado y R-cuadrado ajustado disminuye considerablemente en relación al mismo modelo entrenado en el dataset sin outliers

tidy(anova(modelo_ds_original_ln_base))

Evaluación comparativa del modelo base con el nuevo dataset y el original

Agrupamos los modelos para poder compararlos

# armamos lista con todos los modelos
models_leverage <- list(modelo_ds_original_ln_base = modelo_ds_original_ln_base, modelo_ln_base = modelo_ln_base)
# calculamos las variables resumen
map_df(models_leverage, tidy, .id = "model")

lista_predicciones_training_leverage = map(.x = models_leverage, .f = augment) 

Utilizamos la función glance sobre todos los modelos para comparar la variabilidad de los mismos

df_evaluacion_train_leverage = map_df(models_leverage, glance, .id = "model") %>%
  # ordenamos por R2 ajustado
  arrange(desc(adj.r.squared))
df_evaluacion_train_leverage

Calculamos el RMSE para todos los modelos a evaluar

map_dfr(.x = lista_predicciones_training_leverage, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Calculamos el MAE para todos los modelos a evaluar

map_dfr(.x = lista_predicciones_training_leverage, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Realizando la comparación de las metricas R-cuadrado ajustado, RMSE y MAE, observamos que:

• En la evaluación del nuevo dataset con outliers, el R-cuadrado ajustado disminuyó significativamente por lo que nuestro modelo representa una variabilidad de los datos menor (modelo_ln_base:0.354 - modelo_ds_original_ln_base:0.109)
• El error cuadrático medio aumenta de forma considerable (modelo_ln_base:9.91 - modelo_ds_original_ln_base:15.7)
• El error absoluto MAE, también crece de forma cosiderable (modelo_ln_base:7.46 - modelo_ds_original_ln_base:9.23)

Concluimos con lo observado que nuestro modelo es sensible a outliers, y ante la precencia de estos, nuestras métricas de performance del modelo empeoran drásticamente

Creación modelo linear robusto

Ante la apreciación de como empeora nuestro modelo ante la existencia de outliers en nuestro dataset, procedemos a realizar un entrenamiento con un modelo lineal robusto el cual esperaremos que brinde mejores resultados debido a la ponderación que realizan sobre la influencia de los casos atípicos

# ajustamos modelo lineal múltiple
modelo_ds_original_ln_robust_base <- rlm(peso ~ edad + genero + altura + dias_actividad_fisica_semanal + consumo_diario_alcohol, data = encuesta_salud_train_original_base)
# # Resumen del modelo

tidy(modelo_ds_original_ln_base)

tidy(modelo_ds_original_ln_robust_base)

Comparando coeficientes, podemos ver como los relacionados al modelo robusto, no sufren grandes modificaciones en relación al modelo entrenado con el dataset sin outliers, mientras que en el modelo lineal común entrenado con el dataset con outliers, los coeficientes cambian considerablemente debido a estos datos.

options("scipen"=1)
tidy_ln_ds_original_robust <- tidy(modelo_ds_original_ln_robust_base, conf.int = TRUE)
tidy_ln_ds_original_robust %>% dplyr::select(term, statistic, conf.low, conf.high)

Evaluación comparativa del modelo común y robusto

models_robust <- list(modelo_ln_5_robust = modelo_ds_original_ln_robust_base, modelo_ln_5_base = modelo_ds_original_ln_base)

robust_predictions_over_training = map(.x = models_robust, .f = augment) 
map_dfr(.x = robust_predictions_over_training, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

map_dfr(.x = robust_predictions_over_training, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

• Observamos que el RMSE no mejora con el modelo robusto (16.0)
• Notamos una mejora importante en el error absoluto MAE, enter el modelo lineal (9.23) y robusto (8.76)

Evaluamos la performance de los modelos sobre el conjunto de testing

# Aplicamos la función augment a los 4 modelos con el set de testing
lista_predicciones_testing_robust = map(.x = models_robust, .f = augment, newdata = encuesta_salud_test) 
# Obtenemos el RMSE para los 4 modelos
map_dfr(.x = lista_predicciones_testing_robust, .f = rmse, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

map_dfr(.x = lista_predicciones_testing_robust, .f = mae, truth = peso, estimate = .fitted, .id="modelo") %>% arrange(.estimate)

Evaluando los modelos sobre el conjunto de pruebas

• El error cuadrático medio disminuye con el modelo robusto (10.3) en comparacion con el lineal (10.6) para el dataset con outliers, pero no de forma significativa, no es una métrica robusta, sensible a outliers
• El error absoluto si disminuye con el modelo robusto (7.53) en comparación con el lineal (8.15) de forma significativa para el dataset con outliers, dando un valor similar a la medición del modelo sin la presencia de los valores atípicos

Conclusión, ante la presencia de una cantidad significativa de outliers, utilizar un modelo lineal robusto nos sirvio para poder tener un modelo que ante la presencia de estos valores atípicos se sigue comportando de una manera similar al modelo que fue entrenado sin outliers, mientras que el modelo lineal no robusto, siendo entrenado en un dataset con outliers, cambia significativamente sus coeficientes. Hace sentido poner foco sobre la métria MAE ya que es más robusto cuando los datos tienen outliers o datos atípicos y es la mejor opción a usar en esos casos

Consideraciones

Todos los fragmentos de códigos utilizados en el trabajo práctico fueron orientado a partir de notebooks propuestos por la cátedra En cuanto al modelo robusto donde utilice la librería MAAS, algunas fuentes de ejemplo fueron link1, link2